Parametric Resonance in Bose-Einstein Condensates 
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We demonstrate parametric resonance in Bose-Einstein condensates (BECs) with attractive two- 
body interaction in a harmonic trap under parametric excitation by periodic modulation of the 
s-wave scattering length. We obtain nonlinear equations of motion for the widths of the condensate 
using a Gaussian variational ansatz for the Gross-Pitaevskii condensate wave function. We conduct 
both linear and nonlinear stability analyses for the equations of motion and find qualitative agree- 
ment, thus concluding that the stability of two equilibrium widths of a BEC might be inverted by 
parametric excitation. 
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The phenomenon of parametric resonance, when a sys- 
tem is parametrically excited and oscillates at one of its 
resonant frequencies, is ubiquitous in physics: the phe- 
nomenon is found from simple classical systems such as 
the child's swing and the vertically driven pendulum pQ, 
to the Paul ion trap [5] and aspects of some inflation mod- 
els of the universe [3]. Parametrically excited systems 
are usually nonlinear: even the parametrically driven 
pendulum is governed by the linear Mathieu equation 
only for small oscillations about its equilibrium positions. 
Nonetheless, valuable qualitative insight into these sys- 
tems may be gained by investigating their behavior in 
the vicinity of equilibrium points. 

In the realm of ultracold quantum gases, initial stud- 
ies of parametric resonance have featured Faraday pat- 
terns 4-6 , Kelvin waves of quantized vortex lines [7], 
self-trapped condensates [8j [9], bright and vortex soli- 
tons [101 111) , and self-damping at zero temperature |12j . 
Recently, it was shown for a shaken optical lattice that a 
periodic driving can even induce a quantum phase tran- 
sition from a Mott insulator to a superfluid [H], 
paving the way for new techniques to engineer exotic 
phases (151 116] . Other novel experimental techniques 
[T71 [TB] to excite a Bose-Einstein condensate (BEC) of 
7 Li in the vicinity of a broad Feshbach resonance [H] by 
harmonic modulation of the s-wave scattering length, in 
contrast to the usual method of excitation by modulation 
of the trapping potential [20H23] , have inspired investi- 
gations of parametric resonance and other phenomena 
in Refs. [TTJ [251127] . Conspicuously absent from these 
investigations is a study of the simplest case of para- 
metric resonance in Bose-Einstein condensates: a sys- 
tematic stability analysis for the parametrically excited 
three-dimensional BEC in a harmonic trap. Therefore, 
we conduct in this letter a lucid proof-of-concept study 



of such a parametrically excited BEC by demonstrating 
that within both a linear analytic and a nonlinear nu- 
meric analysis the stability characteristics of its equilib- 
rium configurations can be changed. 

We start with modeling the dynamics of a BEC at zero 
temperature using the mean-field Gross-Pitaevskii (GP) 
Lagrangian 
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The functional derivative of the Lagrangian results in 
the well-known GP equation for the dynamics of the 
mean-field condensate wave function yj = ip(r,t). In ex- 
periments generically a cylindrically-symmetric harmonic 
trapping potential V(r) = mui 2 (p 2 + X 2 z 2 )/2 is used, 
whose elongation is described by the trap anisotropy pa- 
rameter A = u z /u)p. Furthermore, we assume that the s- 
wave scattering length is periodically modulated accord- 
ing to 



a(t) = a + cii sin VLt. 



(2) 



In the following we will investigate how the stability of 
condensate equilibria depends on both the driving am- 
plitude a\ and the driving frequency f2, provided the 
time-averaged s-wave scattering length ao is slightly neg- 
ative. In principle, this could be analyzed by solving the 
underlying GP equation for the condensate wave func- 
tion tp = ip(r,t), which follows from extremizing the La- 
grangian ([I]). However, the thorough numerical analysis 
in Ref. |25| demonstrated convincingly that the dynam- 
ics of the GP equation can be well-approximated within a 
Gaussian variational ansatz for the GP condensate wave 
function [28, 29 1. Even for long evolution times and in 
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the vicinity of resonances, where oscillations of the con- 
densate are quite large, it was possible to reduce the GP 
partial differential equation to a set of ordinary differen- 
tial equations for the variational parameters. 

Therefore, we follow the latter approach and employ 
the Gaussian ansatz 
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with time-dependent variational widths u p , u z , phases 
4> p , 4> z , and normalization M(t) = N 1 / 2 ir 3 / 2 u p 1 u z 1 ^ 2 . In- 
serting the Gaussian ansatz ^ into the GP Lagrangian 
and extremizing with respect to all variational pa- 
rameters, we obtain at first explicit expressions for the 
phases </> PiZ = mu p _ z /(2hu PtZ ). We define the dimen- 
sionless time t — uj p t and scale the variational widths 
by Wp, z = "Mp z/etho, where aho = ^/(mw,,) is the har- 
monic oscillator length. Finally, we write the dimension- 
less driving function p(r) = Po+Pi sin (J7r/w p ) according 
to the definitions po.i = \/'V 7r -^' a o,i/ a ho- The resulting 
dynamics for the widths u p and u z is then determined 
by a pair of coupled nonlinear ordinary differential equa- 
tions: 
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For attractive two-body interactions, there exists a 
critical value of the time-averaged interaction strength 
p" lt {\) < beyond which no equilibria exist in the ab- 
sence of parametric driving. Its dependence on the trap 
anisotropy A must be evaluated numerically, as for ex- 
ample in Ref. [27]- For pQ rlt (A) < p a < 0, there exists a 
pair of equilibrium points for Eqs. Q, one stable and one 
unstable [2SJ [29] , which we denote with u + and u - , re- 
spectively. We remark that the stability of these points in 
the absence of parametric driving is determined by eval- 
uating the frequencies of collective modes for small oscil- 
lations about equilibrium [27] . Whereas the equilibrium 
u + has real frequencies for all modes and, thus, is stable, 
the equilibrium Uo- possesses an imaginary frequency for 
one mode, implying exponential behaviour and thus in- 
stability. 

In view of a linear stability analysis we assume small 
oscillations about an equilibrium, write u p w u p q + 5u p 
and u z ~ u z o + 8u p , and expand the nonlinear terms 
of Eqs. Q to first order in 8u p and Su z . We scale and 
translate time as 2i' + 7r/2 = CIt/uj p , define displacement 
and forcing vectors x(t') and f , 
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and finally we introduce the matrices A and Q corre- 
sponding to constant and periodic coefficients, respec- 
tively: 
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The result is a set of coupled, asymmetric, inhomoge- 
neous Mathieu equations, 

x(i') + (A - 2piQcos2i')x(£') = fcos2t', (7) 

whose solutions determine whether the underlying equi- 
librium is stable or unstable. 

The Mathieu equation, a special case of Hill's differ- 
ential equation |30| . has been studied extensively in the 
literature [31]: approaches to obtaining the equation's 
stability diagram include continued fractions [301 E21 [33J , 
perturbative methods [34H36] . and infinite determinant 
methods 37 . The problem has been treated in great 
detail in Ref. [41] for the study of the Paul trap, the sta- 
bility of which is governed exactly by a set of coupled 
homogeneous Mathieu equations. Of importance to our 
particular problem are Refs. [321 E3], where it was shown 
that for both single and coupled Mathieu equations, an 
inhomogeneous term does not affect the location of sta- 
bility borders to Eq. In the following, we choose the 
concise continued- fraction method based on the approach 
of Refs. [31 [33]. 

The 7r-periodic parametric driving of the Mathieu 
equation permits the application of Floquet theory, the 
essential statement of which is that each of the two fun- 
damental solutions x 12 (£') to Eq. ^ may be written in 
the form |4H 



Xl,2(0 



lint 1 



(8) 



where the 7r-periodic part consists of Fourier components 
b 2n and the exponential part is characterized by the Flo- 
quet exponent /3, which determines the stability of the 
solution. Due to the presence of both signs of the expo- 
nent in a general solution, we require for stability that 
3? [0\ = 0. On the stability borders, solutions to the 
Mathieu equation are m7r-periodic with m£l Thus to 
obtain the stability borders, we set j3 = mi in Eq. (|J). 
By substitution of the Floquet ansatz ([8| into Eq. ([7| , we 
obtain a third-order recurrence relation for the Fourier 
coefficients b 2 „: 



A + (fl + 2in) I b 2n -p 1 Q(b 2 „ +2 + b 2 



0. (9) 
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We define the ladder operators S2„b2 n = i)2n±2 as 

S± = {A + [B + 2i (n + l)] 2 I - PiQSf^y 1 Pl Q, 

(10) 

and by repeated re-substitution of these ladder opera- 
tors into the recursion relation Q for n increasing and 
decreasing from zero, we obtain a tri-diagonal matrix- 
valued continued fraction relating the parameters A, Q, 
and B; 



case 1: Uq- 



case 2: Uq+ 



A + £ 2 I-;p 2 Q 



Q b = 0. (11) 
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In order to obtain a non-trivial solution for b , the 
determinant of the matrix-valued continued fraction in 
Eq. (11) must vanish. The resulting stability for a par- 



ticular value of the modulation frequency fl and the di- 
mensionless driving amplitude pi is then shown in the 
linear stability diagram of Fig. [I] for three characteristic 
values of the trap anisotropy A. Our results indicate that 
for the unstable (stable) equilibrium position, the largest 
region of stability (instability) occurs for a pancake BEC, 
i.e., for A > 1. 

As a special case the stability borders for the isotropic 
condensate are given by a separate continued fraction, 
obtained by an analogous process for a single inhomo- 
geneous Mathieu equation. The resulting stability dia- 
grams are shown in Fig. [2] It allows to draw a direct 
analogy between the isotropic condensate and the para- 
metrically driven pendulum: the pendulum too is de- 
scribed by a single Mathieu equation, however the inho- 
mogeneous term in the equation of motion for the BEC 
corresponds to a direct periodic driving in phase with 
the parametric driving. It was shown in Ref. [42] that 
a periodic inhomogeneity has no effect on the stability 
borders for a single Mathieu equation, so Fig. [2] is simply 
a transformation of the standard Ince-Strutt stability di- 
agram for the parametrically driven pendulum |30) . for 
the relevant experimental parameters of the BEC. 

While the results of linear stability analysis for the cou- 
pled Mathieu system are qualitatively similar to those 
for the single equation, there are a number of notable 
changes between Figs. [T] and [2j First, for the coupled 
Mathieu equations, there exist a set of stability regions 
not attainable by the analytic method used here. These 
are displayed without black borders in Fig. [TJ and corre- 
spond to the so-called "combined resonances" of the sys- 
tem [SHlllI]- These regions are attainable by numerical 
stability analysis of the Mathieu equations [JTJ 0U SS] , 
which was used to generate the colored background re- 
gions of Fig. [T] It is notable in Fig. [I] that these anoma- 
lous regions are not present for A = 1. 

A second and important difference from the single to 
the coupled Mathieu equations is the appearance of a 
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FIG. 1. Linear stability diagrams for the unstable (case 1) 
and stable (case 2) equilibrium positions of a cylindrically- 
symmetric BEC for three values of the trap anisotropy A: (a) 
A = 0.2 (a cigar-shaped BEC), (b) A = 1 (spherical BEC), 
and (c) A = 2.6 (pancake-shaped BEC). White regions corre- 
spond to unstable solutions, darkest shaded regions to stable 
solutions, and lightly shaded regions correspond to marginally 
stable solutions - regions where only one of two available col- 
lective oscillation modes is stable. 



new region, shaded white and issuing from Q/uj p w 10 in 
Fig. [I] for A = 1, case 1. This region is identified with the 
instability of the quadrupole collective mode, which docs 
not appear in a one-dimensional analysis. This result 
implies that a three-dimensional analysis might result in 
further changes to the linear stability diagram of Fig. [I] 
for A = 1, case 1. 

As the underlying equations of motion Q are inher- 
ently nonlinear, a linear analytic stability analysis alone 
is not sufficient to investigate the phenomenon of para- 
metric resonance. Therefore, we have also performed a 
detailed numerical stability analysis by integrating the 
equations of motion Q over a long time-period using a 
Runge-Kutta-Verner 8(9) order algorithm, incrementing 
through pairs (pi , f2) and recording divergent solutions 
to obtain the corresponding stability diagram. The cor- 
responding results are shown in Fig. [3] for the same three 
values of the trap anisotropy A as in Fig. [T] 
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FIG. 2. Linear stability diagrams for (a) unstable and (b) 
stable equilibria of an isotropic BEC. Shaded and white re- 
gions indicate stable and unstable solutions, respectively. The 
presence of a stable region in (a) indicates that an originally 
unstable equilibrium might be stabilized by parametric exci- 
tation. 



The results for the originally stable equilibrium Uo+ 
show both qualitative and even quantitative similarity 
to the linear stability analysis of Fig. [I] i.e., a similar 
tonguelike structure of unstable regions issuing from cer- 
tain points on the vertical axis. The originally unstable 
equilibrium u _ also shows qualitative similarity to the 
linear analysis of Fig. [TJ however the region of stabil- 
ity begins only for much larger modulation frequency f2 
and dimensionless driving amplitude p\. These results 
are both reasonable, as the linear stability analysis is 
only valid for small oscillations - corresponding to large 
(pi, fi) for equilibrium Uo+ and small (pi, f2) for equilib- 
rium uo- • Furthermore, in contrast to Fig. [T] we find 
in the nonlinear stability diagram of Fig. [3] that stability 
is more easily achieved for a cigar-shaped BEC, i.e., for 
A < 1. 

Further comparison of the linear and nonlinear stabil- 
ity diagrams of Figs. [T] and [3] shows the possibility of 
both simultaneous stability of the equilibrium positions, 
and even the possibility of a complete reversal of the 
stability characteristics. In the latter case, the smaller 
equilibrium position would become the only stable width 
of the condensate, which should be experimentally ob- 
servable. A final observation, applicable to the origi- 
nally unstable equilibrium in both linear and nonlinear 
cases, is the existence of a minimum driving amplitude 
pmin necessar y to stabilize the condensate. The value 
p™ 111 = u q (5uq — 1) w 0.17 po is exactly attainable for the 
linear analysis of the isotropic condensate, and in Fig.[T]is 
approximately independent of A in the considered range 
[0.2,2.6]. For the nonlinear analysis, p™ m ps 1.2 po is also 
approximately independent of A. This feature will have 
implications for an experiment, as in conjunction with 
the width of the Feshbach resonance, it dictates the min- 
imum modulation of the applied magnetic field necessary 
to stabilize the condensate. 
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FIG. 3. Nonlinear stability diagrams for the unstable (left) 
and stable (right) equilibria of a cylindrically-symmetric BEC, 
for three values of the trap anisotropy A: (a) A = 0.2 (a cigar- 
shaped BEC), (b) A = 1 (spherical BEC), and (c) A = 2.6 
(pancake BEC). 



We note that any stability diagram depends on the par- 
ticular choice for the time-averaged dimensionless inter- 
action strength pg. The concrete results in Figs. [T}|3] were 
obtained for the particular value po = 0.9 pQ nt (A). As po 
approaches pg rlt (A), we generically observe a growth of 
the stable (unstable) regions in case 1 (2). 

Finally, we conclude that our proof-of-concept investi- 
gation has unambiguously shown that the phenomenon 
of parametric resonance should be experimentally observ- 
able also for Bose-Einstein condensation. However, due 
to the intrinsic nonlinear nature of the underlying GP 
mean-field theory, a linear analysis, like in the Paul trap, 
is not sufficient to quantitatively study the stability di- 
agram. Thus, in order to achieve a destabilization (sta- 
bilization) of a stable (unstable) BEC equilibrium in an 
experiment, a corresponding numerical nonlinear analy- 
sis is indispensable. Regardless, a linear stability analysis 
provides an intuitive and qualitative understanding of the 
physics of parametric resonance in BECs. 

In the present letter we have focused our attention 
upon a periodic modulation of the s-wave scattering 
length around a slightly negative value, which restricts 
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the number of particles in a BEC to the order of a few 
thousand [46] 07]. However, the phenomenon of para- 
metric resonance might be more important for dipolar 
BECs, where in addition to a repulsive short-range and 
isotropic interaction, also a long-range and anisotropic 
dipolar interaction between atomic magnetic or molecu- 
lar dipoles is present. Provided that the dipolar inter- 
action is smaller than the contact interaction, a stable 
dipolar BEC does exist. But a larger dipolar strength 
leads to mutual existence of both a stable and an unsta- 
ble dipolar BEC jJS] whose stability might be changed via 
a periodic modulation of the harmonic trap frequencies 
or the s-wave scattering length. In that context it might 
also be of interest to estimate how quantum fluctuations, 
which are non-negligible for a larger dipolar interaction 
strength [35], change the stability diagram. The case 
for dipolar Fermi gases is probably even more interesting 
from the point of view of parametric resonance, as for 
any dipolar strength a stable equilibrium coexists with 
an unstable one |50j . Thus, parametric resonance may 
offer a simple efficient approach for realizing equilibria of 
dipolar quantum gases whose properties have so far not 
yet been explored. 
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